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DERIVATION OF A THREE DIMENSIONAL NUMERICAL WATER QUALITY 
MODEL FOR ESTUARY AND CONTINENTAL SHELF APPLICATION 


By Malcolm Spaulding 
University of Rhode Island* 


SUMMARY 


A derivation is given for a three dimensional mass transport equation 
which is appropriate for numerical modeling of estuary and Continental Shelf 
water quality variations for Both the time dependent and steady state cases, 

A stable and accurate finite difference approximation to the derived equation 
is presented and a solution scheme for the resulting equations outlined. 
Preliminary results are obtained using the model for extremely simple problems 
which have analytical solutions. The results indicate that the numerical 
model as presented will provide a fruitful scheme to study water quality 
problems in coastal waters for both steady state and time dependent cases. 


INTRODUCTION 


The ability to quantitatively assess the influence of waste discharges on 
the water quality of receiving waters to ensure their proper management has 
progressed rather rapidly in recent years. Development of mathematical or 
numerical water quality models for the one dimensional (l,2,3,^) and two 
dimensional (5,6,7,8,9,10) cases have been shown to provide a basis for 
prediction of various water disposal alternatives and their effect on water 
quality. Summaries of the status of research in this area are provided in 
state-of-the-art reports cited in References (ll) and (12). 

In each of these models, however, an assumption has been made that one 
or several of the spatial dimensions is not significant for the particular 
area under study. For instance, the one-dimensional approximations normally 
assume that the vertical and cross stream structure of a given water quality 
parameter is of secondary importance while the longitudinal or main axis of 
flow directions is of primary interest. Two-dimensional models character- 
istically eliminate either the vertical or lateral structure while maintain- 
ing the other two. This integral or averaging approach has normally been 
taken, since as the number of spatial dimensions increases, the computational 
difficulties, as well' as computer time and storage requirements to solve the 


^Research performed during the NASA-ASEE 1973 Summer Faculty Fellow Program. 
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resulting equations, increases markedly. There are, however, regions both in 
estuaries and on the Continental Shelf, where a useful representation of the 
actual water quality distribution can only be obtained from a three- 
dimensional picture. 

In addition to averaging over the spatial dimensions to eliminate 
spatial fluctuations of water quality indicators, it is also possible to 
average the predictions over some time scale, ranging from the order of a 
portion of a tidal cycle to many days. These models then, generally are 
classified as time-dependent and steady- state, respectively, and each provide 
an important approach to coastal zone pollution studies. 

The present work will derive a three-dimensional mass transport equation 
appropriate for numerical modeling of estuaries and Continental Shelf water 
quality variations for both the time-dependent and steady-state cases, A 
finite difference approximation to the original equation will be shown and a 
solution scheme for the resulting equations outlined. The tidal hydraulics 
or coastal circulation is assumed to be known from field data or hydraulic 
model studies . 


SYMBOLS 


D ,D ,D 
x’ y’ z 


e ,e ,e 
x’ y* z 


g 

H 

J 

L 


coliform bacteria concentration 
source of coliform bacteria 
Chezy coefficient 

dispersion coefficient in x, y, and z directions respectively 
mean sea level depth 

turbulent diffusion coefficients in x, y, and z directions 
respectively 

gravity acceleration 

total depth, mean sea level depth plus tidal height 

discharge rate of B.O.D. 

reaeration coefficient for D.O. 

decay coefficient for coliform bacteria 

B.O.D. decay coefficient 

biological oxygen demand (B.O.D. ) concentration 



dissolved oxygen (D.O.) concentration 

D.O. saturation level 

mass density or concentration vector 

Richardson number 

salinity or source (sink) strength 

source (sink) strength of D.O. 

source (sink) vector 

time 

mean flow speed 

velocities components in x, y, and z direction respectively 
wave height 
wave length 
wave period 

rectangular cartesian coordinates 




To accurately predict the movement of a water quality parameter, 
requires one to quantitatively answer 
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(a) Hov is a given pollutant from a source transported and distributed 
through the receiving water as a function of time? 

(b) How rapidly does the decay or generation by natural processes add 
or subtract from the water quality parameter being employed as an 
indicator? 


The first question is primarily one involving the fluid mechanics of 
mass transport. The processes involved are advective transport of a consti- 
tuent due to the mean tidal velocity and dispersive transport produced by 
turbulent mixing. These variables are related to physical hydrodynamic 
characteristics of the area under consideration and are both time and space 
dependent. The second question is chiefly concerned with chemical and micro- 
biological processes of species generation and decay which are temperatiire, 
time, and concentration dependent. The conceptual separation employed here 
will be more clearly illustrated for specific water quality variables later 
in this report. 

From these elementary considerations, numerical water quality models 
describing both the hydrodynamic transport and decay or generation of a 
specific constituent assume the form of a mass transport equation with a 
specific set of so\irce-sink or reaction matrix terms. The exact form of the 
reaction matrix, as shown in the next section, depends on whether the consti- 
tuent is conservative or nonconservative and whether its magnitude depends 
on other constituent concentrations. Examples of this conceptual approach 
are presented in references (6) and (10). 

In the following paragraphs, the basic mass transport equation will be 
presented and in a later section it will be shown how this can be extended 
to include water quality parameters. 

In a turbulent medium the mass transport equation for a given consti- 
tuent may be written (l3, 1^) as: 


3P, 

3t 


9x 




V 

3z 


_ 3 , . 3 , 


3y 3z 3z 


+ S (2.1) 


where 

- mass density or concentration of substance A 

e ,e , and e - turbulent diffusion coefficients 
x’ y’ z 

u, V, and w - time mean velocity over short sampling times in the x, y, 
and z directions respectively 

S - source or sink of substance A 


k 



In this form, several approximations have already "been introduced. 
Molecular diffusion has been neglected in anticipation that in the environ- 
ments under consideration it is several orders of magnitude smaller than 
the turbulent eddy diffusivity. The diffusion terms are obtained by assuming 
that the turbulent flux terms such as u'p^ can be adequately modeled by the 
product of an eddy diffusion coefficient and the ensemble mean concentration 
gradient. It has also been assumed that no diffusive transports are caused 
by thermal or pressure gradients within the system. The velocities in the 
system represent time averages over short sampling periods much smaller than a 
tidal cycle but larger than the instantaneous variations i.e. one minute 
periods. This procedure then eliminates the stochastic variations in mass 
density. The approximations that have been shown are essentially valid for 
estuary and Continental Shelf waters (l4, 15)* 

A coordinate system (fig, l) can be chosen with z measured vertically 
upward from the bottom or x and y plane where x is the longitudinal 
axis and y is the lateral axis. The water level caused by the ebb and 
flood of the tide, then is a perturbation, above or below the mean sea 
level plane . 

In order to make the mass transport equation more adaptable to the large 
scale depth variations found in estuary and Continental Shelf regions, and 
also to eliminate troublesome numerical boundary conditions at the sea bed 
and free surface it has been found convenient to nondimensionalize the vertical 
axis using the total depth, that is mean sea level depth plus instantaneous 
tidal height. The resulting mass transport equation then becomes: 


3(H p^) 8(H up^) a(H vp^) a(H cop^) 


at 


3x 




an 


- + 1_ fw 4. i_ ffz 

3x 3x 3y 3y 3ri H 3ri 


SH 


( 2 . 2 ) 


where -ir-* and tt” have been assumed small and 
3x 3y 

H - total depth, mean sea level depth plus tidal height (D + ^) 
D - mean sea level depth 
^ - tidal height 
n = z/H 


and the relationship between the real vertical velocity, w, and the 
nondimen sional vertical velocity, O), is given by: 


w - l\\ 


at 


(2.3) 
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where the derivatives are evaluated holding n constant. 

Equation (2,2) then represents a mathematical model for the convection 
and dispersion of any constituent with a generalized source or sink 

combination of strength S, 


Reaction Model 


Water quality models* as previously stated, generally assume the form of 
a mass transport equation as shown in eq. (, 2 . 2 ) with a specific source--sink : 
or reaction matrix term which has been simply noted as S, In the present 
section, the details of this reaction matrix term will be presented for 
several water quality parameters. 

To aid in the discussion, it is noted that water quality parameters can 
be divided into two general classes--conservative and nonconservative 
models. The non-~conservative models can be further subdivided into those 
containing multistage reaction schemes. In order to better understand this 
division, a few examples are presented. 

1. Conservative Case (Salinity) 

3S ^ 3uS ^ 3vS ^ 3wS 3 , 3S. 3 , 3S. 3 / 3Sv _ _ 

3t * ^ * F" 3^' - ^ ^ 'V 

S - salinity or chlorinity (2,i+) 

2. Nonconservative, Single-Stage Reaction Case (Coliform) 

3C ^ 3uC ^ 3vC ^ 3wC 3 . 3C. 3 , 3C. 

F * F" * ^ ~ - F - F ‘V 

C - coliform bacteria concentration 

K, “ decay rate for coliform bacteria 
d 

Cg - source of coliform bacteria 

3. Nonconservative, Two-Stage Reaction Case 


3 / Tf p 

3z z 3z d 


+ C. 


(2.5) 


(Dissolved Oxygen (D.O) and Biochemical Oxygen Demand (B.O.D.) 
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3L . 3uL ^ 3vL ^ 3vL 3 / ^ ^ ^ - 

3t 3x 3y 3z " 3x ^®x 3x' * 3y ^®y 3y^ " 3z ^ z 3z' " 

- Kj^L + J(x,y,z,t) 

12. j. 1 h2 4 . 1x2 4 . 1x2 1_ / 12^ 1_ r 12^ 1_ fp 12\ _ 

9t 8x 3y 9z " 3x ^®x 3x^ " 3y ^®y 3y^ “ 3z ^®z 3z^ “ 

- * ^a(°sat - °> - Sq 

L - B.O.D. concentration 
0 - D.O. concentration 

- B.O.D. decay coefficient 

J - point load of B.O.D. due to outfall 

- reaeration coefficient for D.O, 

OgAT “ saturation level 

Sq - source or sink of D.O. 

From these elementary examples, one can see the essential mass transport 
equation on the left-hand side of the equal sign remaining unchanged in form, 
and the varying source-sink terms on the right-hand side. 

In an attempt to generalize Eq, (2.2) for all water quality parameters, 
we can assume that S, the source-sink term, can be subdivided as: 

S = [K] P + S 


( 2 . 6 ) 


(2.7) 


where S - source or' sink vector 
[K] - reaction matrix 

P - mass density or concentration vector 
then eq. (2.2) may be expressed as 


l£ + li^ + 1x1. + 1x1 

3t 3x 3y 9z 


It gi - |i - fc " J' ■ • 


t ( 2 . 8 ) 


T 



This provides a general equation for all water quality parameters or 
reaction schemes. In its most general form, this approach allows for decay 
or hirth rates dependent on concentration levels of any individual or group 
of constituents, as well as multistage reaction mechanisms. 

As a simple illustration of how this generalization encompasses the 
previous examples, let us look at the reaction matrix for salinity. 


P = [P^] where - salinity 

K = [0] 

S = [0] 


and for the multistage D,.0. - B.O.D, system 


P^ - D.O. concentration 
Pg - B.O.D. concentration 


K 


->'a 

0 -K^ 


s = 


VSAT " Sq 


where the variables have been previously defined. Pulse loads of waste B.O.D. 
defined by J are handled as pulse loads to the sink-source vector $ at 
appropriate spatial locations. 

The generalization of all water quality equations to this conceptual 
form provides a simplification of numerical models since it provides a 
relatively easily handled programming technique to deal with a great variety 
of water quality problems. 
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Dispersion in Multi-dimensional Computations 

In our derivation of. the mass transport equation, it was assumed that the 
normal turbulent mass transport terms such as u'Pa could be adequately 
represented by a diffusion coefficient times the mean concentration gradient. 
It then remains to determine either analytical or empirical expressions that 
can be employed to represent these terms. 

In anticipation of using a finite difference approach for the solution 
of our mass transport equation for which a finite or discrete space is 
assumed to have fixed properties, we must ask ourselves how this approach 
could effect any estimates that are made for the diffusion coefficient in 
question. Physically the discretization process averages the velocities, 
mass densities, and any other variables over the region of our finite space 
increments. Therefore, in any formulation that is made for diffusion 
coefficients we must remember this averaging process has been performed. 

The work of Fischer (l6) provides a straightforward explanation of the 
mechanisms involved in the diffusion of constituents in real environments and 
will provide the basis for determining the coefficients in light of the finite 
difference approximation. In general, various parts of natural water 
environments taken perpendicular to the mean flow show differences in velo- 
city, Due to these variations, portions of a field of pollutant constituent 
will move more or less rapidly than the mean flow, hence dispersing the 
pollutant in the direction of the flow. This process then causes cross 
sectional variations in the mass concentration and leads to a cross sectional 
turbulent diffusion which tends to transfer constituents from parts with 
higher concentration to those with lower concentration. This explanation is 
valid for the contributions from both the vertical and lateral variations 
in the mean flow. Reviewing this description, we have included not only the 
turbulent diffusion but also the effect of shear in the velocity profiles if 
a finite approximation for a spatial grid is used. Under these circumstances 
the "diffusion” coefficients are called dispersion coefficients and will be 
noted in the remainder of this report by D^, D^, and D^ for the x, y, and 
z coordinate directions. 

Elder (iT) based on Taylor’s (l8) concept has determined expressions 
for the longitudinal and lateral dispersion coefficients based on the mean 
velocity, depth, and bottom roughness for simple one-dimensional steady flow 
which should be useful in our problem. Their expressions are; 

D^ = 5.93 Du g^^^ (2,9) 


and 


- 1/2 -1 
D = .23 Du C/ 
y ° z 


( 2 . 10 ) 
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where 


u - mean uniform flow speed 

g - gravitational acceleration 

C - Chezy coefficient 
z 

D - water depth 

Expressions of this type, with a correction of an additive constant for 
wind effects, have been used rather successfully by Leenderste (T) in two- 
dimensional, vertically-averaged, water quality modeling and should provide 
reasonable results in the present work. 

Using an approach based on work by Pritchard et al. (19» 20, 21) in the 
James River estuary, an approximate form of a vertical dispersion coefficient- 
founded on a mixing length theory — which includes both stratification and 
surface mixing due to the wind has been proposed. The formulation is given 
as : 


D 

z 




(D - Z)^ 

Z 



(1 + + 


Z(D - Z) TO 
“P DOT 


2itZ 

WL 


X (1 + SpR^)"^ 


( 2 . 11 ) 


Where 

Tipj^pjOtp adjustable coefficients dependent on particular region of 
interest 

U - mean flow speed 
D - mean sea level depth 

Z - distance measured vertically downward from mean sea level 

- 2 

- Richardson number [g/p 3p/9Z/(.T ] 

WH - wave height 
WL - wave length 
WT - wave period 
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Although this formiilation has only been shown to be strictly applicable 
to the James River area xmder essentially steady-state conditions, it provides 
a good first estimate and could— through the variable constants— probably 
be adapted to give reasonable results for o\ir area of concern. 

As presented, the values of dispersion obtained from various works in 
the literature show that if a comparison of vertical to longitudinal or 
lateral dispersion coefficients is made the former is always several orders 
of magnitude smaller than the latter two. This fact complicates the modeling 
effort since dispersion becomes anisotropic. However, using arguments from 
the work of Holley (22), it can be shown that for most coastal marine environ- 
ments the movement or transport of a particular constituent is determined 
largely by the advection. rather than dispersive transport process in the 
mean sea level plane, and, hence, lateral or longitudinal dispersion is a 
secondary effect. (This is not true in areas immediately surrounding an 
outfall.) Therefore, one can assume an isotropic approximation in each 
direction, and ignore anisotropic effects. 

As a general approach, one finds that as the ability to accurately 
predict the three-dimensional time dependent velocity field increases the 
importance of the dispersion terms decreases. Hence for any particular model 
development , tests should be performed to ascertain the sensitivity to 
changes in dispersion coefficient. This will be mainly dependent on how 
accurate — the size of spatial grids and type of averaging of the tidal 
hydraulic eq.uations— — the velocity or circulation is determined. Therefore, 
the more accurate the circulation picture is known, the less we need to be 
concerned with the dispersive transport mechanism. 


COMPUTATIOKAL MODEL 


Derivation of Finite Difference Annroximation 

Time-Dependent Model . - The mass transport equation as developed in 
equation (2.2) permits a rather large number of finite difference approxi- 
mations. For each possible approach, an analysis of the convergence and 
stability characteristics has to be performed such that ass\irance is gained 
that the difference approximations will actually represent the solution to 
the proposed equation. Following earlier computational work in the field 
of mass transport modeling, notably the work of Leenderste, (6), an extension 
of the A.D.I. (Alternating Direction Implicit) technique was considered as 
the most versatile of the methods since it has unconditional stability 
but only tridiagonal equations to be solved for each direction in space. 

The computational model presented here for the mass transport equation 
is based on the methods developed by Douglas (23), Douglas and Gunn (2U) 
and others ( 25 , 26 ). 

As a first step, a space staggered grid system is chosen to locate the 
discrete values of the variables, (fig. 2. ) It is to be noted that this 
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choice of grid formulation provides centrally-located spatial derivatives 
for the linear terms. In addition, the mass density previously defined as 
which is a point function will be replaced by P indicating that an 
averaging over a physical volume corresponding to the spatial grid size has 
occurred. The following notation for defining discrete values for the 
finite difference approximations will be used. 


X = m Ax 
y = h Ay 
n = n An 
t = A At 


(3.1) 


where m,k,n,£. - integers 

Ax, Ay, An - spatial grid sizes in the x, y, and n coordinate 
directions 


At - temporal 

Therefore, the mass density of a given constituent can be expressed in finite 
difference form as: 


P(x,y,njt) = P(mAx,kAy,nAn,5'At) 


= P 


m,k,n 


(3.2) 


and also defining the difference operators as: 


6 

X m,k,n 

6 P^ 
y m,k,n 


<5 P 

n m,k,n 


m+l/2,k,n m-l/2,k,n 

p'^ _ 

m,k+l/2,n m,k-l/2,n 

P^ - P^ 

m,k,n+l/2 m,k,n-l/2 


(3.3) 


Now employing the approach which Douglas (23) originally used on the 
heat conduction problem, equation (2.2) can be approximated in finite 
difference form by the three following equations: 
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- I « i 6 (UHP)'^ + i 5 (D H5 

2 X 2 X 2 X X X ' 

+ I S (D H6 (P))^ - 6 (VHP)*' - 6 (ojHP)* + 6 (D H5 (P))* 

y ri jjj 

- ^ ^ tSH)* = (z.k) 

advancing the solution from time level Jl to time level Jl + 1/3 and 

- i 6 (UHP)*'^*^^ - i « (UHP)* + i 6 (D H6 (P))***^^ 

2 X 2 X 2 X X X 

+ ^ 6 (D h6 (P))^ - ~ 6 ^ 6 (VHP)^ 

2 X X X ' 2 y 2 y 

“ 6 (D H6 ^ 6 (D H6 (P))^ - 6 (o)HP)^ 

2 y y y 2 y y y r\ 

. yp„. . 

advancing the solution from time level il + l/3 to time level il + 2/3 and 

- ^ 6 (UHP)^'^^'^^ “ o ^ (UHP)^ + ^ (D H6 

2x 2x 2x xx 

+ ^ 6 (D H6 (P))^ - 6 6 (VHP)^ 

2 X X X 2 y 2 y 

.|6y(DyH6y(P))*"2/3,|yc^H6^(p))4 

- 1 ^ i Vr 

, 1 6 MP))* = (3.6) 

2 n H n At 

advancing the solution from time level S. + 2/3 to time level £ + 1. 

An algebraic simplification of the above equations can be obtained by 
subtracting (3.^) from (3.5) and (3.5) from (3.6), respectively. After 
rearrangement the resulting equations become 
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= 5 (UHP)*' - 6 (D H6 (P))®' + 25 (VHP)^ 

X XXX y 

* 25^(^)^ - 26y(DyH6y(P))^ - 25^ (/ 5^(P))^ 


- (HP)®' + 2HS®''^^^^ . 

5 (D H5 (P))®''^^''^ - 5 (VHP - |r (HP)®*^^^ = 

y y y y 

5 (D H5 (P))® - 5„(VHP)® - |- (HP)®*®^^ . 

y y y y 

- 5^(0®)®^® t 5^(^ 5^(P))®^®-|^(HP)®*® = 

-6^(uHP)® + 5^(/-5^(P))®-f^(HP)®^2/3 


(3.7) 


(3.8) 


(3.9) 


Equations (3-7), (3.8), and (3.9) now form a finite difference ^ 
approximation to the original mass transport eq.uation of order 0((Ax) + 

(At)2) and can he solved for any particular rectangular grid system by 
solving the associated tridiagonal linear equations once for each equation. 
References 6 and 10 provide greater details of this approach and show how 
both open and closed boundaries can be handled. 


Steady-State Model .- Since there are many situations when the steady- 
state solution of the mass transport equation (elliptic equation) is of 
particular interest in an area under study, it would be desirable to determine 
necessary modifications to the existing numerical procedure to handle this 
situation. Following the work of Douglas (23) and Wachpress (27) the time 
step increment At in Eqs. (3»7)» (3*8) and (3.9) can be replaced by a 
positive number iteration parameter or sequence of iteration parameters and 
by iteration with these parameters a steady state solution obtained. 

The problem then becomes, after convergence of the solution is assumed, 
the determination of a sequence of iteration parameters which when applied in 
some cyclic pattern will cause the rate of convergence to be maximized. 

Since the literature (23, 2U, 27) provides only an indication of possible 
iteration parameters for a simple heat diffusion problem with constant 
dispersion, an optimum sequence of parameters is not available for the 
general mass transport equation and normally has to be determined through 
numerical experiments. Indications have been made in the work of Aziz and 
Heliums (28) of a possible set of iteration parameters but are not directly 
applicable to this case. 
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Extension of Numerical Model to Include Multistage Water Quality Parameters. - ■ 


Using the reaction matrix approach, as previously outlined, the computational 
schane can he extended to include multistage reactions appropriate to simulate 
various water quality reaction systems by altering the generalized source- sink 
term (HS)^+1/2 as noted in Equation (3.7 )• An appropriate numerical 
approximation then becomes: 

Ul/2 ^ t+1/2 ^ Ul/2 (3.10) 

1 iJ J 1 


where 

i - specific element of the mass density vector P 

.1 - maximum number of constituents 

max 

K - reaction matrix (reaction coefficients) 

S. - sources or sinks of mass density i 
1 

H - total coastal water depth 

Note that P now becomes a mass concentration or density vector containing 
Jmax constituents and each of these must be solved at each complete time 
step. Careful examination of the terms in Equation (3.10) shows that the mass 
density of a given constituent must be evaluated at time level i + 1/2, £,+1/2 

however when solving Equation (3.7) containing this term information for P 
is not available. Douglas (23) has shown that if P is evaluated at the 
lower known time level A-, then the numerical solution has only 0((At) + 

(Ax)^) accuracy. The original accuracy can be regained by estimating a 
value of P^+l/2 ly using Equations (3.7), (3.8) and (3.9) and then going 
back and readvancing the solution in time with the new best estimate for 
pi+1/2^ This procedure, however, essentially requires a doubling of the 
computational time to simulate a fixed physical time, and, therefore, should 
be employed with a consideration of the trade-off between computational 
time and accuracy. 


EFFECTS OF COMPUTATIONAL MODEL APPROXIMATIONS 


Stability and Convergence 

Douglas (23) has shown that for the simple time dependent and steady 
state heat conduction equation the finite difference approximation scheme 
presented yields both a stable and convergent solution with an accuracy of 
0((At)2 + (Ax)^). Indications from other work (28) also indicates that at 
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least in numerical experiments, these properties of convergence and stability 
are still valid when the lower order advection terms are included, as in the 
mass transport equation. Theoretical proof of this case, as far as the 
author knows, has not been determined at this time. 


Mass Conservation 


Another extremely important aspect of any numerical approach to the 
mass transport equation is that mass should be conserved both over the entire 
field under consideration as well as on the near-field or grid point scale. 
Preliminary results with the model as outlined show that for at least the 
simple case of flow in channels of variable depths that this condition is 
adequately met (under 1 percent cumulative mass loss or gain over several 
days s imul at ion ) . 


RESULTS 


To date, the results have been only of a preliminary nature in which 
results for extremely simple problems have been obtained. The time dependent 
model has been run for simple channel cases of both constant and varying 
depths and has shown mass conservation within 0.1 percent over several days 
simulation. The case of simple dispersion — (no advection) has been performed 
for a typical estuary geometry with insignificant mass losses or gains. 
Attempts have also been made to link a two-dimensional vertically averaged 
tidal model to the present mass transport model, but have met with only 
limited success. Although the overall mass conservation appears to be 
within acceptable bounds, local perturbations of mass density near sharply 
varying geometric features cause local variations of the order of 0-4 
percent of the mean field concentration. 

The steady-state model has been rim using a sequence of iteration 
parameters in a cyclic manner for a simple one-dimensional constant velocity, 
constant dispersion, uniform channel flow. The sequence of parameters have 
been chosen based on numerical experiments and is by no means the optimum, 
but provides reasonable convergence rates. 

Figure 3 shows a comparison between the analytical and numerical solu- 
tion to the simple case stated above for several combinations of velocity 
u, dispersion D , and numerical grid spacing Ax. The results were ob- 
tained by 20 cycles of the iteration sequence 


5 X 10 
2,5 X 10“' 

1.25 X lO' 
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6.25 X 10 

3.125 X 10^ 


over a grid field (32 x 12 x 7) in 5 minutes on a GDC 6600 computer with an 
absolute error of 10”^. 


CONCLUSIONS 


A derivation has been given for a three-dimensional mass transport 
equation which is appropriate for numerical modeling of estuary and 
Continental Shelf water quality variations for both the time-dependent and 
steady-state cases. A stable and accurate finite difference approximation 
to the derived equation was presented and a solution scheme for the resulting 
equations outlined. Preliminary results were obtained using the model for 
extremely simple problems which have analytical solutions. These results 
indicate that the numerical model, as presented, will provide a fruitful 
scheme to study water quality problems in coastal waters for both steady- 
state and time-dependent cases. 
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D - Estuary mean sea level depth 
^ - Tidal height 

u.v.w - Velocities in the x,y, and n directions respectively 


P - Mass density 

D ,D , D - Dispersion coefficients for the x,y, and n directions 
^ y 2 respectively 


Figure 2.- Space staggered grid system for 3-D mass transport finite 
difference approximation. 
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C/C (Ratio of concentration to boundary concentration) 


Analytic solution 



Figure 3 - Comparison of numerical steady-state with analytical solution 

for 1-D constant velocity and dispersion. Uniform channe* flow 
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